clear variables; close all; clc;
%% The variables used here are constant over IPER, but change depending setup configuration, stored in storageTank.xsl/SETUP.
% Setup vars will be used in mf?adapt
% Result vars will be retreived from mf?analyze

% Todo:
% - put wells in new sheet, so they can be seen as setupvariables too

saveOutput  = 1;
basename    = 'storageTank';
save name.mat basename;

%%  Load Period data
[perdata.pernams,perdata.pervals]   = getPeriods(basename);
[perdata.perlen]            = perdata.pervals(:,strmatchi('PERLEN', perdata.pernams));
[perdata.milkSign]          = perdata.pervals(:,strmatchi('QMilk',  perdata.pernams));    % Sign of the RECH injection is positive
[perdata.nper]              = length(perdata.perlen);
[perdata.time]              = cumsum(perdata.perlen);
% save perdata pernams pervals perlen milkSign nper time

%%  load setup variables and make setupstructure
% Load values from xsl
[SETUPparnames,SETUPparvals] = getExcelData('storageTank.xls','SETUP','H');
ANALYZEPRETTY	= round(SETUPparvals(:, strmatchi('ANALYZEPRETTY',  SETUPparnames)));
SERIES          = round(SETUPparvals(:, strmatchi('SERIES',         SETUPparnames)));
REUSECONC       = round(SETUPparvals(:, strmatchi('REUSECONC',      SETUPparnames)));
XWALL           = SETUPparvals(:,       strmatchi('SETUP_XWALL',SETUPparnames));
ZWALL           = SETUPparvals(:,       strmatchi('SETUP_ZWALL',SETUPparnames));
MULTIMILK       = SETUPparvals(:,       strmatchi('MULTIMILK',  SETUPparnames));
POREFF          = SETUPparvals(:,       strmatchi('POREFF',     SETUPparnames));
CONDUC          = SETUPparvals(:,       strmatchi('CONDUC',     SETUPparnames));
NSetup          = length(SERIES);

for i=1:length(SERIES); % create setup struct
    setup(i).basename       = basename;
    setup(i).perdata        = perdata;
    
    setup(i).setupVars.SERIES         = SERIES(i);
    setup(i).setupVars.XWALL          = XWALL(i);
    setup(i).setupVars.ZWALL          = ZWALL(i);
    setup(i).setupVars.MULTIMILK      = MULTIMILK(i);
    setup(i).setupVars.ANALYZEPRETTY  = ANALYZEPRETTY(i);
    setup(i).setupVars.POREFF         = POREFF(i);
    setup(i).setupVars.CONDUC         = CONDUC(i);
    setup(i).setupVars.REUSECONC      = REUSECONC(i);
    
    if i~=length(SERIES);
        setup(i).setupVars.REUSENEXTCONC = REUSECONC(i+1);
    else
        setup(i).setupVars.REUSENEXTCONC = 0;
    end
end

%%  Make gridObj for each setup
for i = 1:length(SERIES)
    MW  = 100;	% [cm] Width of model
    MH  = -17;	% [cm]
    D   = 10;   % [cm] thickness of model; dont forget to set well in well-sheet to .5*D
    dx  = 1;	% [cm] grid cell width
    dy  = 5;	% [cm] grid cell length = thickness of model; The model has become a 2D model, Hele Shaw
    dz  = 1;	% [cm] grid cell size vertical ! Cant be less then one?
    
    % horizontal plate injectionBound; impermeable (glass); same for every setup. in cm
    setup(i).grExtra.xhInjB   = [0 3];      setup(i).grExtra.yhInjB   = [0 D];      setup(i).grExtra.zhInjB   = [-3 -2];
    % vertical plate injectionBound; impermeable (glass); same for every setup. in cm
    setup(i).grExtra.xvInjB   = [2 3];      setup(i).grExtra.yvInjB   = [0 D];      setup(i).grExtra.zvInjB   = [-3 -1];
    
    xGr = [0:dx:MW MW 0:8      setup(1).grExtra.xhInjB setup(1).grExtra.xvInjB 21:25];	% [cm] grid, making sure that MW is always included irrespective of dx !
    yGr = [0:dy:D                   setup(1).grExtra.yhInjB setup(1).grExtra.yvInjB];   % [cm] grid
    zGr = [0:-dz:MH MH              setup(1).grExtra.zhInjB setup(1).grExtra.zvInjB];	% [cm] grid, making sure that 0 is always include irrespective of dz ! 0 is groundlevel
    
    gr  = gridObj(xGr,yGr,zGr);  % create gridobj
    setup(i).gr         = gr;
    
    %%  Create extra index variables within grid
    setup(i).grExtra.ixWall     = round(interp1(gr.xm,1:length(gr.xm),setup(i).setupVars.XWALL));     % Convert to IBOUND, so nodegrid is used
    setup(i).grExtra.izWall     = round(interp1(gr.zm,1:length(gr.zm),setup(i).setupVars.ZWALL));     % Convert to IBOUND, so nodegrid is used
    setup(i).grExtra.ixRech     = round(interp1(gr.xm,1:length(gr.xm),setup(i).setupVars.XWALL));     % Ceil/floor consturtion is used, so RECH is always on the leftside of the wall
    
    % horzontal plate injectionBound
    setup(i).grExtra.ixhInjB	= find(gr.xm >= setup(1).grExtra.xhInjB(1) & gr.xm <= setup(1).grExtra.xhInjB(2));
    setup(i).grExtra.iyhInjB	= find(gr.ym >= setup(1).grExtra.yhInjB(1) & gr.ym <= setup(1).grExtra.yhInjB(2));
    setup(i).grExtra.izhInjB	= find(gr.zm >= setup(1).grExtra.zhInjB(1) & gr.zm <= setup(1).grExtra.zhInjB(2));
    [setup(i).grExtra.iXHInjB,setup(i).grExtra.iYHInjB,setup(i).grExtra.iZHInjB] = meshgrid(setup(1).grExtra.ixhInjB,setup(1).grExtra.iyhInjB,setup(1).grExtra.izhInjB);
    % vertical plate injectionBound
    setup(i).grExtra.ixvInjB    = find(gr.xm >= setup(1).grExtra.xvInjB(1) & gr.xm <= setup(1).grExtra.xvInjB(2));
    setup(i).grExtra.iyvInjB    = find(gr.ym >= setup(1).grExtra.yvInjB(1) & gr.ym <= setup(1).grExtra.yvInjB(2));
    setup(i).grExtra.izvInjB    = find(gr.zm >= setup(1).grExtra.zvInjB(1) & gr.zm <= setup(1).grExtra.zvInjB(2));
    [setup(i).grExtra.iXVInjB,setup(i).grExtra.iYVInjB,setup(i).grExtra.iZVInjB] = meshgrid(setup(1).grExtra.ixvInjB,setup(1).grExtra.iyvInjB,setup(1).grExtra.izvInjB);
end

%%  save and clear other variables
% set outputfolder
if saveOutput
    outputPath=sprintf('%s\\output\\%s\\',pwd,datestr(now,30));
    mkdir(outputPath);
    
    
    for i = 1:length(SERIES),
        setup(i).outputPath = outputPath;
        setup(i).saveOutput = saveOutput;
        setup(i).setupNumb  = i;                % for now, setup(n) == setupNumb, better that way
    end
    save(strcat(outputPath,'setup'),'setup');
end

save setup setup

% setup = matfile(setup,'Writable',isWritable); future option, better for
% large setup structures
%%  Loop through the series
for i=1:length(setup);
    n=i;
    save curNum n
    fprintf('going to mf_setup setup number: %f\n',n)
    
    mf_build;
    mf_setup;
    fillResults; %     interpretate concentration files
end

clear variables